

cap log close	
local gruppo = 1
	
log  using  .../A_8_`gruppo'_2_26Sep2020.smcl,replace

version 13	
use .../sample_estimation_`gruppo'_2_13Sep2020.dta, clear
		
*** 1 ***	
	cap program drop prog2_lnf
	set type double					

	program define prog2_lnf
*       version 8.2

	args todo b lnf

	tempvar ISE IU ISELF IULF logitp2 logitp3 logitp4 b1 b2 b3 b4   a1 a2 a3 a4 logitp5 logitp6 logitp7 logitp8  lnfie last lE1 lE2 lE3   lE1LF lE2LF lE2LF ///
	        lU1  lU2 lU3 lU1LF lU2LF lU3LF l1 ll1 l2 ll2 l3 ll3 l4 ll4 l5 ll5 l6 ll6 l7 ll7 l8 ll8  mlSE1 mlSE2  mlSE3  mlE1 mlE2 mlE3 mlSELF1  mlSELF2 mlSELF3  mlELF1 mlELF2 mlELF3 ///
	        h1SE h2SE h3SE  sum1SE sum2SE  sum3SE  ST1SE ST2SE ST3SE lE3LF h1SELF h2SELF  h3SELF  sum1SELF sum2SELF  sum3SELF  ST1SELF ST2SELF ST3SELF ///
		h1U h2U  h3U  sum1U sum2U  sum3U  ST1U ST2U ST3U h1ULF h2ULF   h3ULF  sum1ULF sum2ULF  sum3ULF  ST1ULF ST2ULF ST3ULF lastSE lastU lastSELF lastULF lnfi
		
	tempname p2 p3 p4 p1
	
        *defining equations
    mleval `ISE'     = `b', eq(1)
	mleval `IU'      = `b', eq(2)
    mleval `ISELF'   = `b', eq(3)
	mleval `IULF'    = `b', eq(4)
	mleval `a1'       = `b', scalar eq(5)
	mleval `a2'       = `b', scalar eq(6)
	mleval `b1'       = `b', scalar eq(7)
	mleval `b2'       = `b', scalar eq(8)
	mleval `logitp2' = `b', scalar  eq(9) //5



	

	quietly {
			* to make sure probabilities constrained between 0 and 1
			* for identification, p11==1-p21-p22 and m11==0
			scalar `p2' = exp(`logitp2')/(1+exp(`logitp2'))
			scalar `p1' = 1             /(1+exp(`logitp2'))

			***** for SELFEMPLOYED spells
			
			* hazard functions for the three types
			by $S_E_id $spell: gen double `h1SE' = 1-exp(-exp(`ISE' + `a1'))	if $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `h2SE' = 1-exp(-exp(`ISE' ))			if $SE_E ==1 & $LF == 0
			
                        * first part of survival function 
			by $S_E_id $spell: gen double `sum1SE' = sum(exp(`ISE' + `a1'))		if $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `sum2SE' = sum(exp(`ISE' ))			if $SE_E ==1 & $LF == 0
		
			
			* cumulative product of survival functions
			by $S_E_id $spell: gen double `ST1SE' = exp(-`sum1SE'[_N])			if _n==_N & $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `ST2SE' = exp(-`sum2SE'[_N])			if _n==_N & $SE_E ==1 & $LF == 0
			

			* identify last obs of the spell
			by $S_E_id $spell: gen byte `lastSE' = (_n==_N) if $SE_E ==1 & $LF == 0
			replace `lastSE'=0 if `lastSE'==.

			* likelihood functions of the spells for each type
			by $S_E_id $spell: gen double `lE1' =  `ST1SE' * ((`h1SE'/(1-`h1SE'))^$event ) if _n==_N & $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `lE2' =  `ST2SE' * ((`h2SE'/(1-`h2SE'))^$event ) if _n==_N & $SE_E ==1 & $LF == 0
			
			
			replace `lE1'=0 if `lE1'==.
			replace `lE2'=0 if `lE2'==.
			

			******** for EMPLOYED 
			by $S_E_id $spell: gen double `h1U' = 1-exp(-exp(`IU' + `b1' )) if $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `h2U' = 1-exp(-exp(`IU' )) if $SE_E ==0 & $LF == 0
			
			
			by $S_E_id $spell: gen double `sum1U' = sum(exp(`IU' + `b1' )) if $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `sum2U' = sum(exp(`IU' )) if $SE_E ==0 & $LF == 0
			

			by $S_E_id $spell: gen double `ST1U' = exp(-`sum1U'[_N]) if _n==_N & $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `ST2U' = exp(-`sum2U'[_N]) if _n==_N & $SE_E ==0 & $LF == 0
			
			
			by $S_E_id $spell: gen byte `lastU' = (_n==_N) if $SE_E ==0 & $LF == 0
			replace `lastU'=0 if `lastU'==.

			by $S_E_id $spell: gen double `lU1' =  `ST1U' * ((`h1U'/(1-`h1U'))^$event ) if _n==_N & $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `lU2' =  `ST2U' * ((`h2U'/(1-`h2U'))^$event ) if _n==_N & $SE_E ==0 & $LF == 0
                        
			               
			replace `lU1'=0 if `lU1'==.
			replace `lU2'=0 if `lU2'==.
			
			
			***** for SELFEMPLOYED left-censored spells
			
			* hazard functions for the three types
			by $S_E_id $spell: gen double `h1SELF' = 1-exp(-exp(`ISELF' + `a2'))   if $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `h2SELF' = 1-exp(-exp(`ISELF'))   if $SE_E ==1 & $LF == 1
			
			
                        * first part of survival function
			by $S_E_id $spell: gen double `sum1SELF' = sum(exp(`ISELF' + `a2'))  if $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `sum2SELF' = sum(exp(`ISELF' ))    if $SE_E ==1 & $LF == 1
			
			
			* cumulative product of survival functions
			by $S_E_id $spell: gen double `ST1SELF' = exp(-`sum1SELF'[_N]) if _n==_N & $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `ST2SELF' = exp(-`sum2SELF'[_N]) if _n==_N & $SE_E ==1 & $LF == 1
			
			
			* identify last obs of the spell
			by $S_E_id $spell: gen byte `lastSELF' = (_n==_N) if $SE_E ==1 & $LF == 1
			replace `lastSELF'=0 if `lastSELF'==.

			* likelihood functions of the spells for each type
			by $S_E_id $spell: gen double `lE1LF' =  `ST1SELF' * ((`h1SELF'/(1-`h1SELF'))^$event ) if _n==_N & $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `lE2LF' =  `ST2SELF' * ((`h2SELF'/(1-`h2SELF'))^$event ) if _n==_N & $SE_E ==1 & $LF == 1
			
           
			replace `lE1LF'=0 if `lE1LF'==.
			replace `lE2LF'=0 if `lE2LF'==.
			
			

			******** for EMPLOYED left-censored spells
			by $S_E_id $spell: gen double `h1ULF' = 1-exp(-exp(`IULF' + `b2')) if $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `h2ULF' = 1-exp(-exp(`IULF' )) if $SE_E ==0 & $LF == 1
			

			by $S_E_id $spell: gen double `sum1ULF' = sum(exp(`IULF' + `b2')) if $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `sum2ULF' = sum(exp(`IULF' )) if $SE_E ==0 & $LF == 1
			
			
			by $S_E_id $spell: gen double `ST1ULF' = exp(-`sum1ULF'[_N]) if _n==_N & $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `ST2ULF' = exp(-`sum2ULF'[_N]) if _n==_N & $SE_E ==0 & $LF == 1
			
			
			by $S_E_id $spell: gen byte `lastULF' = (_n==_N) if $SE_E ==0 & $LF == 1
			replace `lastULF'=0 if `lastULF'==.

			by $S_E_id $spell: gen double `lU1LF' =  `ST1ULF' * ((`h1ULF'/(1-`h1ULF'))^$event ) if _n==_N & $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `lU2LF' =  `ST2ULF' * ((`h2ULF'/(1-`h2ULF'))^$event ) if _n==_N & $SE_E ==0 & $LF == 1
			
                                       
			replace `lU1LF'=0 if `lU1LF'==.
			replace `lU2LF'=0 if `lU2LF'==.
			
	    
				// individual log likelyhood for each type
			by $S_E_id: gen double `l1'=`lE1' +`lU1'+`lU1LF'+`lE1LF'
			by $S_E_id: gen double `ll1'=`p1'*exp(sum(log(`l1')))
			by $S_E_id: gen double `l2'=`lU2'+`lE2' +`lU2LF'+`lE2LF'
			by $S_E_id: gen double `ll2'=`p2'*exp(sum(log(`l2')))




			by $S_E_id: gen byte `last' = (_n==_N) 
	
			* individual contribution to the log likelihood 
			gen double `lnfie' = cond( !`last'  , 0, $weight * ln(( `ll1') + ( `ll2')   ))
				
					
			mlsum `lnf' = `lnfie'

	}
end

	
	
	

set type double		


tempvar mysamp 
		tempname b b01 b02 b1 b2 b03 b04 b3 b4 lnf V  a1 a2 a3  a4 b4 masshE masslE masshSE masslSE ///
		       logitp2 logitp3 logitp4  masslSELF masslELF  masshSELF masshELF

  *now put all the variables togetherlocal
 *   se_netincdiff se_convexity se_female_1 se_ageb se_edu_2-se_edu_3  se_yr_2-se_yr_4 se_reg_2-se_reg_5
 *    lnt_se t_se t2_se 
 *    e_netincdiff e_convexity e_female_1 e_ageb e_edu_2-e_edu_3  e_yr_2-e_yr_4 e_reg_2-e_reg_5
 *  lnt_e  t_e t2_e
 *  self_netincdiff self_convexity  self_female_1 self_ageb self_edu_2-self_edu_3 self_yr_2-self_yr_4 self_reg_2-self_reg_5
 * t_self t2_self lnt_self
 * elf_netincdiff elf_convexity elf_female_1 elf_ageb elf_edu_2-elf_edu_3  elf_yr_2-elf_yr_4 elf_reg_2-elf_reg_5
 * t_elf t2_elf lnt_elf
global varse " se_netincdiff se_convexity se_female_1 se_ageb se_edu_2-se_edu_3  se_yr_2-se_yr_4 se_reg_2-se_reg_5"
global vare  " e_netincdiff e_convexity e_female_1 e_ageb e_edu_2-e_edu_3  e_yr_2-e_yr_4 e_reg_2-e_reg_5"
global varself " self_netincdiff self_convexity  self_female_1 self_ageb self_edu_2-self_edu_3 self_yr_2-self_yr_4 self_reg_2-self_reg_5"
global varelf  "elf_netincdiff elf_convexity elf_female_1 elf_ageb elf_edu_2-elf_edu_3  elf_yr_2-elf_yr_4 elf_reg_2-elf_reg_5"
global dumse "  lnt_se"
global dume "    lnt_e"
global dumself " lnt_self "
global dumelf " lnt_elf"

 
set more off	
// initial values for ll function taken from clogclog		
glm durvar  $varse  $dumse if sedum==1 & lcens_spell ==0 ,  f(b) l(cloglog) 
matrix `b01' = e(b)
matrix coleq `b01' = hazardE
	
glm durvar   $vare  $dume  if sedum==0 & lcens_spell ==0   ,  f(b) l(cloglog) 	
matrix `b02' = e(b)
matrix coleq `b02' = hazardU

glm durvar  $varself  $dumself  if sedum==1 & lcens_spell ==1 ,  f(b) l(cloglog) 
matrix `b03' = e(b)

matrix coleq `b03' = hazardELF
	
glm durvar   $varelf  $dumelf  if sedum==0 & lcens_spell ==1 ,  f(b) l(cloglog) 	
matrix `b04' = e(b)

matrix coleq `b04' = hazardULF

local LL1 = e(ll)
// create matrix with initial betas from cloglog and initial mass point location and probability (given)

matrix `a1' = (-0.1 )
matrix colnames `a1' = a1:_cons

matrix `a2' = (-1.3)
matrix colnames `a2' = a2:_cons

matrix `a3' = (-3.2)
matrix colnames `a3' = a3:_cons

matrix `b1' = (-0.2 )
matrix colnames `b1' = b1:_cons

matrix `b2' = (-1.2 )
matrix colnames `b2' = b2:_cons
matrix `b3' = (-1.5 )
matrix colnames `b3' = b3:_cons


matrix `logitp2' = (0.5)
matrix colnames `logitp2' = logitp2:_cons
matrix `logitp3' = (-2.1)
matrix colnames `logitp3' = logitp3:_cons
matrix `logitp4' = (-1.1)
matrix colnames `logitp4' = logitp4:_cons


matrix `b1' = `b01',`b02',`b03',`b04',`a1',`a2',`b1',`b2',`logitp2'


matrix list `b1'

// define macro that identify id variable
global S_E_id "id"
// define macro that identify se or E spell
global SE_E "sedum"
// define macro that identy the spell
global spell "a_spell"
// define macro that identify event
global event "durvar" 
 // define macro that identify left_censored spells
global LF "lcens_spell" 
 
// define weight
global weight "pweight"
 
 
sort id a_spell year
 

//ml model, d0, no gradient or Hessian provided
ml model d0 prog2_lnf (hazardE: durvar = $varse  $dumse ) ///
                       (hazardU: durvar =  $vare  $dume  ) ///
		       (hazardELF: durvar = $varself $dumself   ) ///
		       (hazardULF: durvar = $varelf $dumelf  ) ///
					    (a1: ) (a2: ) (b1: ) (b2: )    ///
					    (logitp2: )   ///
                                , waldtest(0) /*init(`b1')*/   


ml init hazardE:se_netincdiff    		=  -.4853675 
ml init hazardE:se_convexity    		=   .0157164 
ml init hazardE:se_female_1    			=   .0408632 
ml init hazardE:se_ageb    				=  -.0110562 
ml init hazardE:se_edu_2    			=  -.0372353 
ml init hazardE:se_edu_3    			=   .2093057 
ml init hazardE:se_yr_2    				=  -.0609402 
ml init hazardE:se_yr_3    				=   .0535824 
ml init hazardE:se_yr_4    				=  -.3071054 
ml init hazardE:se_reg_2    			=   .0717283 
ml init hazardE:se_reg_3    			=   .0143267 
ml init hazardE:se_reg_4    			=   .0235397 
ml init hazardE:se_reg_5    			=   .1055123 
ml init hazardE:lnt_se    				=  -.4999314 
ml init hazardE:_cons    				=  -1.237143 
ml init hazardU:e_netincdiff    		=    1.82182 
ml init hazardU:e_convexity    			=   -.319216 
ml init hazardU:e_female_1    			=   .5740778 
ml init hazardU:e_ageb    				=   .0266474 
ml init hazardU:e_edu_2    				=   .1480045 
ml init hazardU:e_edu_3    				=   .2139789 
ml init hazardU:e_yr_2    				=   .2172805 
ml init hazardU:e_yr_3    				=   .2903141 
ml init hazardU:e_yr_4    				=  -.1955885 
ml init hazardU:e_reg_2    				=  -.2781331 
ml init hazardU:e_reg_3    				=  -.3939244 
ml init hazardU:e_reg_4    				=  -.3966443 
ml init hazardU:e_reg_5    				=  -.4916966 
ml init hazardU:lnt_e    				=  -.4304544 
ml init hazardU:_cons    				=  -2.900852 
ml init hazardELF:self_netincdiff    	=  -.6329412 
ml init hazardELF:self_convexity    	=   .0151313 
ml init hazardELF:self_female_1    		=   .1060201 
ml init hazardELF:self_ageb    			=  -.0386727 
ml init hazardELF:self_edu_2    		=  -.1180198 
ml init hazardELF:self_edu_3    		=   .0401427 
ml init hazardELF:self_yr_2    			=   -.166973 
ml init hazardELF:self_yr_3    			=  -.2632063 
ml init hazardELF:self_yr_4    			=   -.781877 
ml init hazardELF:self_reg_2    		=    .026081 
ml init hazardELF:self_reg_3    		=   .0345907 
ml init hazardELF:self_reg_4    		=   .1257196 
ml init hazardELF:self_reg_5    		=   .0917355 
ml init hazardELF:lnt_self    			=  -.0297469 
ml init hazardELF:_cons    				=  -.6244635 
ml init hazardULF:elf_netincdiff    	=    1.75979 
ml init hazardULF:elf_convexity    		=   -.177054 
ml init hazardULF:elf_female_1    		=   .7178054 
ml init hazardULF:elf_ageb    			=   -.042944 
ml init hazardULF:elf_edu_2    			=   .0823624 
ml init hazardULF:elf_edu_3    			=   .2197798 
ml init hazardULF:elf_yr_2    			=   .1711268 
ml init hazardULF:elf_yr_3    			=   .0798423 
ml init hazardULF:elf_yr_4    			=   -.599842 
ml init hazardULF:elf_reg_2    			=   .4285149 
ml init hazardULF:elf_reg_3    			=   .0814714 
ml init hazardULF:elf_reg_4    			=   .0432518 
ml init hazardULF:elf_reg_5    			=  -.0938249 
ml init hazardULF:lnt_elf    			=  -.2609266 
ml init hazardULF:_cons    				=  -2.330838 
ml init a1:_cons    					=  -.6037221
ml init a2:_cons    					=  -1.342118  
ml init b1:_cons    					=   -2.79554 
ml init b2:_cons    					=  -2.067829 
ml init logitp2:_cons    				=  -1.609524
			
ml maximize, trace  search(off) difficult gradient

eststo m

eststo r1lr: nlcom (p1: 1/(1+exp([logitp2]_cons))) ///
                    (p2:exp([logitp2]_cons)/(1+exp([logitp2]_cons))), post

log close

set more off
